Simulation of Multicellular Tumor Spheroids Growth Dynamics 
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The inverse geometric approach to the modeling of the growth of circular objects revealing required 
features, such as the velocity of the growth and fractal behavior of their contours, is presented. 
It enables to reproduce some of the recent findings in morphometry of tumors with the possible 
implications for cancer research. The technique is based on cellular automata paradigm with the 
transition rules selected by optimization procedure performed by the genetic algorithms. 
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Introduction 



Understanding the fundamental laws driving the tu- 
mor development is one of the biggest challenges of con- 
temporary science Q, 0. Internal dynamics of a tu- 
mor reveals itself in a number of phenomena, one of 
the most obvious ones being the growth. Overtaking 
its control would have profound impact to therapeutic 
strategies. Cancer research has developed during the 
past few decades into a very active scientific field tak- 
ing on the concepts from many scientific areas, e. g., 
statistical physics, applied mathematics, and nonlinear 
science USEE Ell E E El 13 El 13 El From 
a certain point of view, the evolution of tumors can be 
understood as an interplay between the chemical inter- 
actions and geometric limitations mutually conditioning 
each other. Consequently, it is believed that malignancy 
of a tumor can be inferred exceptionally from the geo- 
metric features of its interface with the surrounding it 
tissue [Tiil ITU . Formation of the growing interface is in a 
continuum approximation described by a variety of al- 
ternative growth models, such as Kardar-Parisi-Zhang 
|18| , molecular beam epitaxy (MBE) [Tgj . or Edwards- 
Wilkinson equations. At the same time, the growth 
properties can be classified into universality classes [2fj| . 
each of them showing specific scaling behavior with cor- 
responding critical exponents. As implies scaling analysis 
of the 2-dimensional tumor contours [2lt |22| , the tumor 
growth dynamics belongs to the MBE universality class 
characterized by, (1) a linear growth rate (of the radius), 

(2) the proliferating activity at the outer border, and 

(3) diffusion at the tumor surface. 

In vitro grown tumors usually form 3- (or 2-) dimen- 
sional spherical (or close to) aggregations, called mul- 
ticellular tumor spheroids (MTS) |23j. These provide, 
allowing strictly controlled nutritional and mechanical 
conditions, excellent experimental patterns to test the 



validity of the proposed models of tumor growth |24j . 
These are usually classified into two groups, (1) contin- 
uum, formulated through differential equations, and (2) 
discrete lattice models, typically represented by cellular 
automata |H I2fll2^l , agent-based , and Monte Carlo 
inspired models [29j. 

Here we present the inverse geometric approach to the 
MTS growth simulation, enabling us to evolve an initial 
MTS by required rate as well as desired fractal dimension 
of its contour. The method is based on the cellular au- 
tomata paradigm with the transition rules found by the 
genetic algorithms. 



Simulation and optimization tools 

Cellular automata (CA) |3C| were originally introduced 
by John von Neumann as a possible idealization of bio- 
logical systems. In the simple case they consist of a 2D 
lattice of cells s'ff e {0, 1}, i,j = -(L D /2), -(Ln/2) + 
I,.., Ld/2, where t = 0, . . . ,r, is the time step and Ld 
size of the 2D lattice. During the r time steps they evolve 
obeying the set of local transition rules (CA rules) <r, for- 
mally written 
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which defines the CA rules a as the mapping 
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{0,1}x{0,1}x...x{0,1}^{0,1}. (2) 



Any deterministic CA evolution is represented by the 
corresponding point in 2 9 -dimcnsional binary space en- 
abling, in principle, immense number of 2 possible 
global behaviors, predestining CA to be the efficient sim- 
ulation and modeling tool |3l| . Inherent nonlincarity of 
CA models is, however, a double-edged sword. On the 
one hand, it enables to model a broad variety of behav- 
iors, from trivial to complex, on the other hand it results 
in difficulty with finding the transition rules generating 
the desired global behavior. No well-established universal 
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technique exists to solve the problem, and, despite spo- 
radic promising applications of genetic algorithms (GA) 
to solve the task [^2, HE EH > one typically implements 
CA by ad hoc or microscopically reasoned definition of 
the transition rules. 

Genetic Algorithms ( GA ) |35j | are general-purpose 
search and optimization techniques based on the anal- 
ogy with Darwinian evolution of biological species. In 
this respect, the evolution of a population of individ- 
uals is viewed as a search for an optimum (in general 
sense) genotype. The feedback comes as the interaction 
of the resulting phenotype with environment. Formal- 
izing the basic evolutionary mechanisms, such as muta- 
tions, crossing-over and survival of the fittest, the fun- 
damental theorem of GA was derived {schema theorem) 
which shows that the evolution is actually driven by mul- 
tiplying and combining good (quantified by an appropri- 
ate objective function), correlations of traits (also called 
schemata, or hyperplanes). The remarkable feature re- 
sulting from the schema theorem is the implicit paral- 
lelism stating that by evaluating a (large enough) popu- 
lation of individuals, one, at the same time, obtains the 
information about the quality of many more correlations 
of the traits. 

Bellow we present the application of GA optimization 
to find the CA rules producing the 2D CA evolution by 
required rate as well as fractal behavior of the contour. 



Optimization Problem 

In the below numerical studies two competing hypothe- 
ses of the rate of the desired tumor mass production have 
been distinguished, 

a) broadly accepted exponential growth of the tumor 
mass: 



M« = B + 



tt(R 
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exp(Ci) , 



(3) 



and, 

b ) th e growth of the mass with linearly growing radius 



MW = n(RW + At) 2 



(4) 



where R^ is the initial cluster radius; A,B,C are con- 
stants parameterizing the growth process. The term 
B + [-7r(i? (0) ) 2 — B] in Q was chosen to start from the 
initial cluster size = n(R^) 2 for any B. 

At the beginning, the chain of concentric circles (pat- 
terns) with randomly deformed close-to-circular contours 

pg } £ {0, 1}, i,j = -(Ld/2), (-Ld/2) + 1, . . . , L D /2, for 
t = 0, . . . , r, were generated accordingly to 
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, or 1 drawn with probability 
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1 , for yji 2 + j 2 < flW - 1, 



(5) 



where the increasing tumor radius is taken as 



(6) 



in the case of exponential tumor mass production Eq. , 
and 



R (t) = R (o) + M 



(7) 



in the case of linearly growing radius model Eq. re- 
spectively. 

The optimization task solved by GA was to find 
the CA transition rules a* Eq. @ providing the 



growth from the initial pattern {sS } = {pfj'}, hj — 
— (Ld/2), . . . , (Ld/2) — l,Lx>/2 through the sequence of 
the square lattice configurations {sy }, t = 1, . . . ,r, gen- 
erated accordingly to Eq. {Tjl. with the minimum differ- 
ence from the above "prescribed" patterns {p« } in the 
respective t, quantified by the objective function 
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where 5 is the standard Kronecker delta symbol, wq is the 
weight factor, in our case wq = 2. The above expression 
of the objective function Eq. @ reflects the programming 
issues. The larger overlap of {s\f} with {p|*- } enhances 
the denominator of Eq. (JHJl, the prefactor pij in the term 
PijS PijSij reduces the computational overhead by ignoring 
the calcul S PljSlj for p tj = 0. 

The other requirement to the desired growth relates 
to the geometric properties of the contour/interface. 
Broadly accepted invariant measure expressing the con- 
tour irregularity is the fractal dimension, Dp. Using the 
box-counting method it can be calculated as 



Dv — lim 



\ogN B (e) 



o log(l/e) 



(9) 



where N-^(e) is the minimum number of boxes of size e 
required to cover the contour. Here, it has been deter- 
mined as the slope in the log- log plot of Nb(£) against 
1/e using the standard linear fit. 

To obtain the CA rules generating the cluster with 
the required fractal dimension (er*), Dp, the objective 
function Eq. © has been multiplied by the factor 



f 2 (a) = l + w 1 (D F -D ¥ 



(10) 



where D F is the fractal dimension of the cluster kept after 
the r steps with the CA rules tr; the weight w\ was kept 
1 in all the presented numerical results. 

Finally, the rule-dependent objective function has been 
written 



f(a)=h(a)f 2 (a). 



(11) 



To sum up, the optimum rule a* is the subject of GA 
optimization 



t/(<0 = /(O- 



(12) 
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Results and discussion 

All the CA runs started from the pattern {p^ '} 
Eq. (J5J, with the radius = 5. In all the below 

GA optimizations, all the CA evolutions ran on the 2- 
dimensional box of the size Ld = 300 cells and length of 
the CA evolution r = 100 (300) steps. The GA search has 
been applied to find the set of CA rules, a*, which gives 
minimum objective function values (Eq. (JSJ, or Eq. Q11JI- 
respectively). The size of the population was kept con- 
stant (1000 individui), the probability of bit-flip mutation 
0.001, the crossing over probability 1, and ranking selec- 
tion scheme applied. The length of the optimization was 
3000 generations. 




100 150 200 

CA steps 



Exponential growth vs. linear radius dilemma 

The simplest mathematical models of MTS growth [3(J 
assume exponential increase of the MTS mass during the 
time Eq. @. The above assumption is applied mainly 
because of feasibility of differential and integral calcu- 
lus, nevertheless revealing, hopefully, some of the char- 
acteristics of the real growth. Bru et al. pll 12*2^ have 
shown from the morphometric studies that the mean ra- 
dius of 2D tumors grows linearly. At the same time, they 
have experimentally shown that the cells proliferation is 
located near the outer border of cell colonies. In this 
work the former assumption has been tested and com- 
pared with the alternative exponential growth (Figs. ^ 
|5J). For that reason the GA search has been carried out to 
find the CA transition rules which minimize the criterion 
Eq. (JSJ) with exponentially growing pattern Eq. JSJ during 
r = 100 steps. In the inset of Fig. ^one can see that on 
the interval of the optimization (r = 100) the CA evolu- 
tion can be in principle fitted by the exponential Eq. Q 
(as it was required) as well as by Eq. |0J, corresponding 
to the growth with linearly increasing radius, both with 
small systematic error, which can be possibly hidden in 
the stochasticity of the real biological growth. On the 
other hand, the extrapolation of the fits beyond the in- 
terval of optimization shows obvious divergence of the 
CA mass production from the exponential fit, meanwhile 
its deviation from the regime with linearly growing ra- 
dius stays small (nevertheless systematic). We attribute 
the latter discrepancy to the fact that the growing inter- 
face during the CA evolution is neither smooth, nor of 
zero thickness, which is true also for real tumors growth. 
In Fig. [21 we show the comparison of the CA mass pro- 
duction, using the same CA rules as in Fig. 2] fitted by 
the Eqs. © and (@J, respectively, on the interval much 
longer than the interval of minimization (t = 100). One 
can see that both the fits are still possible, nevertheless 
the exponential fit deviates crucially on the interval of op- 
timization (the inset of Fig. [2Jl ■ The above results show 
better agreement of the CA mass production with the hy- 
pothesis of the MTS growth with linearly growing radius 
Eq. with a slight implication towards experiments - 



FIG. 1: Comparison of the CA mass production using the 
CA transition rules obtained by the minimization of Eq. (0 
during r = 100 steps against the exponential growth Eq. 
(dashed line) and the growth with linearly increasing radius 
Eq. (solid line). The extrapolation of both the fits beyond 
the interval of minimization is displayed. The inset shows 
coincidence of the fits on the interval of the minimization. 




CA steps 

FIG. 2: Comparison of the CA mass production using the 
same CA transition rules as in Fig Here, the CA mass 
production has been fitted by both the fits far beyond the 
interval of optimization. The inset shows the comparison of 
the fits on the interval of optimization. 



the growth of close-to-circular MTS is probably slightly 
faster than proposed by Bru et al. [^I, nevertheless not 
exponential. 



Fractal behavior of the contour 

Figs. |3] and 0] show the efficiency of the above approach 
to simulate the MTS growth by any required rate (Eq. (J5J 
with coming from Eq. Q), and, at the same time, 
desired fractal dimensions of the contour, Dp. Here, two 
GA optimizations have been performed to find the CA 
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rules generating the mass production minimizing both 
the criteria Eq. (JSJ and Eq. I|1U|) in the multiplicative 
form Eq. during r = 300 steps and reaching the 
desired fractal dimensions D F — 1.35 (Fig. |3J), and 1.1 
(Fig. 0J, respectively. The obtained CA rules provide 
the growth that fits well to the required rate as well as 
desired D-p. Fig. [3] shows the final size of the CA clus- 




FIG. 3: The CA mass production using the CA rules obtained 
by GA minimization of the Eq. 11 U with the patterns Eq. JSJ 
substitute by Eq. J7J with i?' ' = 5 and velocity constant 
A — 0.4, and the desired fractal dimension Dp 00 = 1.35 (here 
presented CA run gives Df° = 1.34). The smooth circles 
correspond to Eq. 




FIG. 4: The same as in Fig- 01 requiring the fractal dimension 
Dp 00 = 1.1 (here presented CA run gives Dp 00 = 1.09). 

ters using the CA rules obtained by the GA minimization 
of Eq. (|TT)l requiring the growth accordingly to Eq. Q 
during the r = 100 steps within a broad range of the 
velocity constants, A E< 0.1,0.5 >, as well as the frac- 
tal dimensions, £> F 00 G< 1.0,1.35 >. Fig. shows the 



average Dp 00 generated by the above CA rules for each 
of the respective pairs of the parameters A, Dp. The 
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FIG. 5: Resulting contours after r = 100 steps using the CA 
rules found by the GA for the respective pairs of the desired 
parameters A, Df- The smooth circles correspond to the 
growth accordingly to Eq. |g) for the respective A. The real 
values of Dp 00 that correspond to the above CA clusters after 
r = 100 steps are depicted in Fig. [S] 




FIG. 6: The real fractal dimensions Dp 00 using the CA rules 
found by GA optimization for the respective pairs of A, Df 
(Fig. . Dotted pi ane shows respective desired fractal dimen- 
sions, Df. 

above results demonstrate that a specific fractal behavior 
of the growing interface (Fig. [SJ) out of the broad range 
(Dp G< 1.05, 1.35 >) consistent with the morphometric 
results |2lJ can be intentionally generated by here pre- 
sented GA optimization approach. Moreover, the growth 
rate can still be kept at desired value. Beyond these lim- 
its unwanted artifacts emerge. 
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Scaling behavior 

Bru et al. 0, used the scaling analysis to char- 
acterize the geometric features of the contours of a few 
tens growing tumors and cell colonies. Here we outline 
the scaling behavior of the contour of the growing CA 
cluster resulting from our approach. 

A rough interface is usually characterized by the fluc- 
tuations of the height h(x, t) around its mean value h, 
the global interface width [37| 

W(L,t) = ([h(x,t) -h\ ) , (13) 

the overbar is the average over all x, L is an Euclidean 
size and the brackets mean average over different realiza- 
tions. In general, the growth processes are expected to 
follow the Family- Vicsek ansatz |3g| . 

W(L,t)=lPf(L/S(t)), (14) 

with the asymptotic behavior of the scaling function 

f^) = \ const if u»l, ( 15 ) 

where a is the rougheness exponent and characterizes the 
stationary regime in which the height-height correlation 
length f (i) - t 1 / 2 (z is so called dynamic exponent) has 
reached a value larger than L. The ratio (3 = a/z is 
called the growth exponent and characterizes the short- 
time behavior of the interface. 

To adapt the scaling ansatz to the close-to-circular 
growing CA cluster we identify the constant Euclidean 
size L with the time- dependent mean radius r = r(t) = 

(T,k=i } r k (t))/N c (t), r k (t) being the distance of the k- 
th contour point from the center and Nc(t) the numcr 
of contour points in t. Subsequently, we rewrite Eq. 1|13|) 
into 

W c (t) = ({r k (t)-rf) 1/2 , (16) 

the overbar being the average over all the contour points 
in t and the brackets mean average over different realiza- 
tions of contours reaching the mean radius r in t, with 
the scaling ansatz Eq. 114(1 applied 

W c {t)=t?f{r/m)- (17) 

Numerical investigation of the Wc (t) (Fig. [7J| reveals the 
region with the power law behavior. We identify the 
respective slope in the log-log plot with the growth ex- 
ponent = 0.723 (fitted in the region 150 < t < 700). 

To draw more complete scenario of the growing CA 
cluster scaling behavior, the deeper investigation of site- 
site correlation functions in both radial and poloidal di- 
rections is needed. This type of studies will follow. 
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FIG. 7: Scaling behavior of Wc(t). The region with power law 
behavior can be identified, corresponding to the slope 0.723. 



Conclusion 

In the paper we have presented the approach to the 
modeling of multicellular tumor spheroids by required 
growth rate and fractal dimension. The technique is 
based on the combination of the CA modeling with the 
transition rules searched by the GA. Here demonstrated 
results show the feasibility of the approach in this specific 
case. Based on the similarity of the geometric properties 
of the CA evolution and the tumor contours (such as lo- 
cality of the interaction/communication, deformed con- 
tour and nonzero thickness of the proliferating layer) we 
have reasoned that the MTS mass production is slightly 
faster than corresponding to linearly growing radius |22| . 
At the same time our results imply that the often used 
Gompertz curve of the tumor mass progression comes as 
a higher level phenomena related to the nutrition, space 
restrictions, etc. We believe that our approach could 
be implemented as the backbone into the more sophisti- 
cated models of tumor growth encompassing the above 
higher-level mechanisms. Computationally efficient on 
the fly scaling analysis during the CA evolution would 
enable to bias the GA optimization towards the MTS 
growth with desired scaling properties of the contour; 
nevertheless its efficient realization is the subject of our 
ongoing research. Successful implementation of scaling 
analysis into the optimization process could significantly 
contribute to the discussion on the scaling behavior of 
the real tumors [2 I l39t l40| . 

In our opinion the above presented optimization ap- 
proach to the modeling of growing clusters by required 
rate and surface properties can find applications in many 
different fields, such as molecular science, surface design 
or bioinformatics. 
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